% Theoretical EXIT chart ploting

clear all

EbN0_db_vect=1/2;

R=1/3;
mes_size=10000;

g=[1 1 1; 1 0 1];
% g=[0 1 0 0 1 1;0 1 1 1 1 1];
% generator polynomials parameters
[n,K] = size(g);
% number of memories
m = K-1;

muA_vect=[0.1:0.5:2 3 6 9 20];

IA=zeros(length(EbN0_db_vect),length(muA_vect));
IE=zeros(length(EbN0_db_vect),length(muA_vect));
for j=1:length(EbN0_db_vect)
    EbN0_db=EbN0_db_vect(j);
    
    [random, pos] = sort(rand(1,mes_size+m));
    Ave=5;
    for k=1:Ave
        x= round(rand(1, mes_size));
        TurboEncoderOutputs = TurboEncode(g, x, pos);

        encodedBits= zeros(1,3*length(TurboEncoderOutputs(1,:)));
        for i=1:mes_size
            encodedBits(3*i-2:3*i) = (TurboEncoderOutputs(:,i)');
        end
        encodedSymb = 2*encodedBits - 1;

        EbN0=10^(EbN0_db/10);

        sigmaB_2= 1/(2*R*EbN0);
        Lc=2/sigmaB_2;
        bruit= randn(1,length(encodedSymb)).*sqrt(sigmaB_2);
        noisySymb= encodedSymb + bruit;
        noisyDecIn = [ noisySymb(1:3:end);
            noisySymb(2:3:end);
            noisySymb(3:3:end)];
        i=0;

        for i= 1:length(muA_vect)
            muA=muA_vect(i);
            
            clear pAp;
            clear pAm;
            clear pEp;
            clear pEm;
            clear LaMinus;
            clear LaPlus;
            clear extrMinus;
            clear extrPlus;
            
            sigma_A=sqrt(2*muA);
            nA=randn(1,mes_size+m)*sigma_A;
            La= muA*(2*TurboEncoderOutputs(1,:)-1) + nA;

            LaPlus= La(find(TurboEncoderOutputs(1,:)>0.5));
            LaMinus= La(find(TurboEncoderOutputs(1,:)<0.5));

            pAp=hist(LaPlus,min(La):1:max(La));
            pAp=pAp/sum(pAp);% pAp(x) = pA(x|X=1)

            pAm=hist(LaMinus,min(La):1:max(La));
            pAm=pAm/sum(pAm);% pAm(x) = pA(x|X=-1)

            IA(j,i)=IA(j,i)+1/2*(sum(pAp(find(pAp>0)).*log2(2*pAp(find(pAp>0))./...
                (pAp(find(pAp>0))+pAm(find(pAp>0)))))...
                +sum(pAm(find(pAm>0)).*log2(2*pAm(find(pAm>0))./...
                (pAp(find(pAm>0))+pAm(find(pAm>0))))));

            [LMAP extrinsic] = TurboDecode1(g, noisyDecIn(:,1:end),Lc, mes_size+m, La');

            extrPlus= extrinsic(find(TurboEncoderOutputs(1,:)>0.5));
            extrMinus= extrinsic(find(TurboEncoderOutputs(1,:)<0.5));

            pEp=hist(extrPlus,min(extrinsic):1:max(extrinsic));
            pEp=pEp/sum(pEp);% pEp(x) = pE(x|X=1)

            pEm=hist(extrMinus,min(extrinsic):1:max(extrinsic));
            pEm=pEm/sum(pEm);% pEm(x) = pE(x|X=-1)

            IE(j,i)=IE(j,i)+1/2*(sum(pEp(find(pEp>0)).*log2(2*pEp(find(pEp>0))./...
                                        (pEp(find(pEp>0))+pEm(find(pEp>0)))))...
                                +sum(pEm(find(pEm>0)).*log2(2*pEm(find(pEm>0))./...
                                        (pEp(find(pEm>0))+pEm(find(pEm>0))))));
        end
    end
end
IA=IA/Ave;
IE=IE/Ave;

figure;hold on
xx = 0:.05:1;
for i=1:length(EbN0_db_vect)
    yy = spline(IA(i,:),IE(i,:),xx);
    plot(IA(i,:),IE(i,:),'o',xx,yy,IE(i,:),IA(i,:),'r*',yy,xx,'r')
end

legend('','IE1(IA1)','','IA2(IE2)')
xlabel('IA1,IE2')
ylabel('IE1,IA2')





